# Indicators
allIndicators <- read_csv(here::here("Processed_Indicators/allIndicators.csv"))
allIndicators$Season <- factor(allIndicators$Season, levels= c("spring", "summer", "fall", "winter", "all"))
allIndicators$stat_area <- factor(allIndicators$stat_area, levels= c("511", "512", "513", "511-513"))
indicators <- unique(allIndicators$stat_area)
plot_fun <- function(x, indicator){
x %>%
filter(Indicator == indicator) %>%
ggplot() +
geom_line(aes(Year, Value, col = stat_area)) +
theme_bw() +
facet_wrap(~Season)
}
sumLms_fun <- function(indicator, season, stat_areas){
df <- allIndicators %>%
filter(Indicator == indicator,
Season == season,
stat_area == stat_areas)
lm1 <- lm(Value ~ Year, data = df)
return(lm1)
}
acf_plot <- function(indicator, season, stat_areas){
df <- allIndicators %>%
filter(Indicator == indicator,
Season == season,
stat_area == stat_areas) %>%
select(Year, Value)
acf(df)
}
pscore_fun <- function(indicator, season, stat_areas){
df <- allIndicators %>%
filter(Indicator == indicator,
Season == season,
stat_area == stat_areas)
lm1 <- lm(Value ~ Year, data = df)
pscore <- segmented::pscore.test(lm1)
return(pscore)
}
davies_fun <- function(indicator, season, stat_areas){
df <- allIndicators %>%
filter(Indicator == indicator,
Season == season,
stat_area == stat_areas)
lm1 <- lm(Value ~ Year, data = df)
davies <- segmented::davies.test(lm1)
return(davies)
}
slope_bp_fun <- function(indicator, season, stat_areas, Npsi){
df <- allIndicators %>%
filter(Indicator == indicator,
Season == season,
stat_area == stat_areas)
lm1 <- lm(Value ~ Year, data = df)
segmented::segmented(lm1, seg.Z = ~Year, npsi = Npsi)
}
mean_bp_fun <- function(indicator, season, stat_areas, models){
df <- allIndicators %>%
filter(Indicator == indicator,
Season == season,
stat_area == stat_areas)
lm1 <- lm(Value ~ Year, data = df)
temp_mcp1 <- mcp::mcp(models, data = df, par_x = "Year")
return(temp_mcp1)
}
This report includes the initial analysis of the ecosystem and habitat indicators. Linear regression was used to assess trends over the length of the time series. A trend breakpoint analysis and a mean breakpoint analysis were run on the indicators. A breakpoint found in the slope of the line indicates a change in the trend of the data. A breakpoint in the mean of the date may indicate a regime shift or change of the overall state of the system. These analyses were run seasonally within statistical areas, when available, and as yearly averages across the entire study domain.
The following indicators are used in this report
plot_fun(allIndicators, "oisst")
lms <- sumLms_fun("oisst", "all", "511-513")
plot(lms)
summary(lms)
##
## Call:
## lm(formula = Value ~ Year, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.15193 -0.32238 0.06949 0.35912 1.38051
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -78.585742 9.579605 -8.203 3.80e-12 ***
## Year 0.039408 0.004789 8.230 3.38e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4944 on 78 degrees of freedom
## Multiple R-squared: 0.4648, Adjusted R-squared: 0.4579
## F-statistic: 67.73 on 1 and 78 DF, p-value: 3.376e-12
acf_plot("oisst", "all", "511-513")
The pscore test and davies tests indicate whether a breakpoint in the slope exists. The davies test is more sensitive to data that follow a more sinusoidal curve but less sensitive to linear changes. A significant p-value indicates there is a breakpoint in the data. A non-significant p-value indicates a breakpoint is not present.
pscore_fun("oisst", "all", "511-513")
##
## Score test for one/two changes in the slope
##
## data: formula = Value ~ Year
## breakpoint for variable = Year
## model = gaussian , link = identity , method = lm
## observed value = 3.3271, n.points = 10, p-value = 0.00134
## alternative hypothesis: two.sided (1 breakpoint)
davies_fun("oisst", "all", "511-513")
##
## Davies' test for a change in the slope
##
## data: formula = Value ~ Year , method = lm
## model = gaussian , link = identity
## segmented variable = Year
## 'best' at = 1994, n.points = 8, p-value = 0.0005195
## alternative hypothesis: two.sided
bp <- slope_bp_fun("oisst", "all", "511-513", 1)
summary(bp)
##
## ***Regression Model with Segmented Relationship(s)***
##
## Call:
## segmented.lm(obj = lm1, seg.Z = ~Year, npsi = Npsi)
##
## Estimated Break-Point(s):
## Est. St.Err
## psi1.Year 1992 2.029
##
## Meaningful coefficients of the linear terms:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 86.40086 52.38220 1.649 0.103
## Year -0.04359 0.02637 -1.653 0.102
## U1.Year 0.10377 0.02738 3.790 NA
##
## Residual standard error: 0.4459 on 76 degrees of freedom
## Multiple R-Squared: 0.5757, Adjusted R-squared: 0.559
##
## Convergence attained in 1 iter. (rel. change 3.4373e-08)
plot(bp)
df <- allIndicators %>%
filter(Indicator == "oisst",
Season == "all",
stat_area == "511-513")
points(x = df$Year, y = df$Value)
mean_cp <- mean_bp_fun("oisst", "all", "511-513", list(Value~1, 1~1))
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 80
## Unobserved stochastic nodes: 4
## Total graph size: 580
##
## Initializing model
summary(mean_cp)
## Family: gaussian(link = 'identity')
## Iterations: 9000 from 3 chains.
## Segments:
## 1: Value ~ 1
## 2: Value ~ 1 ~ 1
##
## Population-level parameters:
## name mean lower upper Rhat n.eff
## cp_1 2009.354 2008.53 2.0e+03 1.2 1198
## int_1 -0.073 -0.19 4.1e-02 1.0 4238
## int_2 1.089 0.91 1.3e+00 1.0 2597
## sigma_1 0.427 0.36 5.0e-01 1.0 4734
plot(mean_cp)
plot_fun(allIndicators, "fvcom_sst")
lms <- sumLms_fun("fvcom_sst", "all", "511-513")
plot(lms)
summary(lms)
##
## Call:
## lm(formula = Value ~ Year, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.35607 -0.35722 -0.03152 0.43646 1.60338
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -64.956743 10.689314 -6.077 3.94e-08 ***
## Year 0.032429 0.005345 6.068 4.10e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5726 on 80 degrees of freedom
## Multiple R-squared: 0.3152, Adjusted R-squared: 0.3066
## F-statistic: 36.82 on 1 and 80 DF, p-value: 4.097e-08
acf_plot("fvcom_sst", "all", "511-513")
The pscore test and davies tests indicate whether a breakpoint in the slope exists. The davies test is more sensitive to data that follow a more sinusoidal curve but less sensitive to linear changes. A significant p-value indicates there is a breakpoint in the data. A non-significant p-value indicates a breakpoint is not present.
pscore_fun("fvcom_sst", "all", "511-513")
##
## Score test for one/two changes in the slope
##
## data: formula = Value ~ Year
## breakpoint for variable = Year
## model = gaussian , link = identity , method = lm
## observed value = 3.7145, n.points = 10, p-value = 0.0003754
## alternative hypothesis: two.sided (1 breakpoint)
davies_fun("fvcom_sst", "all", "511-513")
##
## Davies' test for a change in the slope
##
## data: formula = Value ~ Year , method = lm
## model = gaussian , link = identity
## segmented variable = Year
## 'best' at = 1988.9, n.points = 8, p-value = 0.0003583
## alternative hypothesis: two.sided
bp <- slope_bp_fun("fvcom_sst", "all", "511-513", 1)
summary(bp)
##
## ***Regression Model with Segmented Relationship(s)***
##
## Call:
## segmented.lm(obj = lm1, seg.Z = ~Year, npsi = Npsi)
##
## Estimated Break-Point(s):
## Est. St.Err
## psi1.Year 1989.496 1.977
##
## Meaningful coefficients of the linear terms:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 163.65082 80.13934 2.042 0.0445 *
## Year -0.08267 0.04038 -2.047 0.0440 *
## U1.Year 0.13529 0.04105 3.296 NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.5187 on 78 degrees of freedom
## Multiple R-Squared: 0.4521, Adjusted R-squared: 0.431
##
## Convergence attained in 2 iter. (rel. change 0)
plot(bp)
df <- allIndicators %>%
filter(Indicator == "fvcom_sst",
Season == "all",
stat_area == "511-513")
points(x = df$Year, y = df$Value)
mean_cp <- mean_bp_fun("fvcom_sst", "all", "511-513", list(Value~1, 1~1))
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 82
## Unobserved stochastic nodes: 4
## Total graph size: 594
##
## Initializing model
summary(mean_cp)
## Family: gaussian(link = 'identity')
## Iterations: 9000 from 3 chains.
## Segments:
## 1: Value ~ 1
## 2: Value ~ 1 ~ 1
##
## Population-level parameters:
## name mean lower upper Rhat n.eff
## cp_1 2009.47 2009.02 2010.00 1 4233
## int_1 -0.43 -0.54 -0.31 1 5679
## int_2 0.79 0.61 0.97 1 5404
## sigma_1 0.43 0.37 0.50 1 5211
plot(mean_cp)
plot_fun(allIndicators, "fvcom_bt")
lms <- sumLms_fun("fvcom_bt", "all", "511-513")
plot(lms)
summary(lms)
##
## Call:
## lm(formula = Value ~ Year, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.95338 -0.35913 0.03703 0.33016 1.27851
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -27.638215 8.590416 -3.217 0.00187 **
## Year 0.013804 0.004295 3.214 0.00189 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.4602 on 80 degrees of freedom
## Multiple R-squared: 0.1144, Adjusted R-squared: 0.1033
## F-statistic: 10.33 on 1 and 80 DF, p-value: 0.001888
acf_plot("fvcom_bt", "all", "511-513")
The pscore test and davies tests indicate whether a breakpoint in the slope exists. The davies test is more sensitive to data that follow a more sinusoidal curve but less sensitive to linear changes. A significant p-value indicates there is a breakpoint in the data. A non-significant p-value indicates a breakpoint is not present.
pscore_fun("fvcom_bt", "all", "511-513")
##
## Score test for one/two changes in the slope
##
## data: formula = Value ~ Year
## breakpoint for variable = Year
## model = gaussian , link = identity , method = lm
## observed value = 2.7247, n.points = 10, p-value = 0.007903
## alternative hypothesis: two.sided (1 breakpoint)
davies_fun("fvcom_bt", "all", "511-513")
##
## Davies' test for a change in the slope
##
## data: formula = Value ~ Year , method = lm
## model = gaussian , link = identity
## segmented variable = Year
## 'best' at = 1993.3, n.points = 8, p-value = 0.0244
## alternative hypothesis: two.sided
bp <- slope_bp_fun("fvcom_bt", "all", "511-513", 1)
summary(bp)
##
## ***Regression Model with Segmented Relationship(s)***
##
## Call:
## segmented.lm(obj = lm1, seg.Z = ~Year, npsi = Npsi)
##
## Estimated Break-Point(s):
## Est. St.Err
## psi1.Year 2004 3.813
##
## Meaningful coefficients of the linear terms:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 10.064436 17.203649 0.585 0.560
## Year -0.005135 0.008636 -0.595 0.554
## U1.Year 0.053291 0.018968 2.810 NA
##
## Residual standard error: 0.4404 on 78 degrees of freedom
## Multiple R-Squared: 0.2093, Adjusted R-squared: 0.1789
##
## Convergence attained in 1 iter. (rel. change -1.3997e-08)
plot(bp)
df <- allIndicators %>%
filter(Indicator == "fvcom_bt",
Season == "all",
stat_area == "511-513")
points(x = df$Year, y = df$Value)
mean_cp <- mean_bp_fun("fvcom_bt", "all", "511-513", list(Value~1, 1~1))
## Compiling model graph
## Resolving undeclared variables
## Allocating nodes
## Graph information:
## Observed stochastic nodes: 82
## Unobserved stochastic nodes: 4
## Total graph size: 594
##
## Initializing model
summary(mean_cp)
## Family: gaussian(link = 'identity')
## Iterations: 9000 from 3 chains.
## Segments:
## 1: Value ~ 1
## 2: Value ~ 1 ~ 1
##
## Population-level parameters:
## name mean lower upper Rhat n.eff
## cp_1 2009.22 2008.00 2010.187 1 3279
## int_1 -0.19 -0.30 -0.085 1 5191
## int_2 0.40 0.22 0.573 1 4080
## sigma_1 0.41 0.35 0.479 1 5006
plot(mean_cp)